
function [y] = Monmer (npoint, hstart, hend, dhg, freq)
    %UNTITLED2 Summary of this class goes here
    %   Detailed explanation goes here

    global filex

    global nlines 
    global hint exper
    global title
    global r1 r2 c1 c2 obox
    global maxx maxy 
    global xy
    global filex
    global nspex ishape nthe nphim theu phiu
    global gx gy gz dx dy dz spinm axas ayas azas niso wgt
    global pop arat nnit npro axbs aybs azbs
    global gxt gyt gzt dxt dyt dzt spinmt axat ayat azat
    global popt aratt nnitt nprot axbt aybt azbt nisot wgtt theut  phiut
    global isaveh ilin iiso
    global flx fly flz flxt flyt flzt bwidemax 

    
    
    iangle=0;
    if(nthe==1 && nphim==1) iangle=1; end
   
    % Need to calculate lineshape further into wings for Lorentians
    if(ishape==1 || ishape==4) 
      out=3.0;
    else
      out=6.0;
    end
    
    runtime=1;
    
    y(npoint)=0.0;
    
    %conversion to cm-1, and account for 1e9 for freq in GHz
    wo=0.03356*freq;
    
    while(runtime==1)
    
        nline=2*spinm+1;
        sp2=spinm*(spinm+1);
        hint(1)=1;
        nhy=1;
        if(nnit>0) nhy=2*nnit+1; end
        if(npro>0) nhy=npro+1; end

        %Now calculate factors that are independent of orientation
        if(ishape==3 || ishape==6) 
          flx2=flx^2;
          fly2=fly^2;
          flz2=flz^2;
        end
        if(ishape==7)
          %parameters related to gaussian distribution for C60
          bz=wo/(beta*gz);
          bx=wo/(beta*gx);
          bmax=mx+2.0*bwidemax;
          bdifmax=(bx-bz)+2.0*bwidemax;
          bdifmxhf=0.5*bdifmax;
        end
        if(spinm~=0.0)
          axa=1.0E-4*axa;
          aya = 1.0E-4*ayas;        %**************************HOW DO YOU USE AYAS BEFORE IT'S INITIALIZED???
          aza = 1.0E-4*azas;
          anon2 = (aya^2-axa^2)*(aya^2-axa^2);
          axya2 = axa^2*aya^2;
        end
        if(nnit~=0 || npro~=0)
          axb=1.0E-4*axbs;
          ayb = 1.0E-4*aybs;
          azb = 1.0E-4*azbs;
        end
        gxy2=gz^2*gy^2;
        gxz2=gx^2*gz^2;
        gyz2=gy^2*gz^2;
        if(iangle==0)
          % half of the steps are in equal increments of theta up to 43 
          % degrees and then rest are in equal steps of cos theta, applies only
          % to powder calculations
          ntheh=nthe/2;
          % calculate increment in theta up to 43 degrees
          step1=.7505/ntheh;
          % calculate increment in cos theta for 43-90 degrees
          step2=.7314/ntheh;
        end

        % Start of theta loop
        % write(*,*) 'Start of theta loop'

        thep=.7505;

        for(it=1: nthe)
            if(iangle==1) 
                the=theu*0.01745;
                costhe=cos(the);
                delth=1;
            elseif(it>ntheh)
                itt=it-ntheh;
                costhe= .7314-(2*itt-1)*step2/2;
                the=acos(costhe);
                delth=the-thep;
                if(itt==1) 
                   delth= 2*(the-.7505); 
                end
                thep=the;
            else 
                the=(2*it-1)*step1/2;
                costhe=cos(the);
                delth=step1;
            end
            csthe2=costhe^2;
            snthe2=1-csthe2;
            sinthe=sqrt(snthe2);
            st2ct2=snthe2*csthe2;
            %Calculate number of phi increments for this value of theta for
            %powder calculations
            if(iangle==1)
                nphi=1;
                step=1.0;
            else
                nphi=nphim*sinthe+1;
                step=1.5708/nphi;
            end
            %start of phi loop
            for(ip=1:nphi)
               if(iangle==1)
                   phi=phiu * .01745;
               else
                   phi=(2*ip-1)*step/2;
                   if(nphi==1) phi=0.0;
                   end
               end
               cosphi=cos(phi);
               sinphi=sin(phi);
               csphi2=cosphi^2;
               snphi2=sinphi^2;

               % Now evaluate terms involving theta and phi

               st2cp2=snthe2*csphi2;
               st2sp2=snthe2*snphi2;
               fgx=gx^2*st2cp2;
               fgy=gy^2*st2sp2;
               fgz=gz^2*csthe2;
               csq= fgx + fgy + fgz;
               g= sqrt(csq);
               if(spinm~=0.0)
                   ak2=(axa^2*fgx + aya^2*fgy + aza^2*fgz)/gsq;
                   ak=sqrt(ak2);
                   if(iangle==1) 
                       nlines=nlines+1 ;
                       akout=1.0e4*ak;
                       %*********************  Write something
                   end
               end
               if(nnit~=0 || npro~=0)
                   akb2 =(axb^2*fgx + ayb^2*fgy + azb^2*fgz)/gsq;
                   akb = sqrt(akb2);
               end

               tpro= (gxy2*snthe2+gyz2*(snphi2+csthe2*csphi2)+gxz2*(csphi2+csthe2*snphi2))/(gsq*g*2);
               if(iangle==1) 
                    tprob = tpro;
               else
                   tprob = tpro*sinthe*delth*step;
               end

               gperp2= gx^2*csphi2+gy^2*snphi2;
               if(spinm~=0)
                   aperp2=(axa^2*gx^2*csphi2*aya^2*gy^2*snphi2)/gperp2;
               end

               d2=(dx^2*fgx+dy^2*fgy+dz^2*fgz)/gsq;
               d=sqrt(d2);    

               % Calculation of orientation dependent fraction lorentzian
               if(ishape==3 || ishape==6) 
                    fl = sqrt((flx2*fgx+fly2*fgy + flz2*fgz)/gsq);
               end

               betg = beta*g;
               betg2 = betg*betg;
               ho=wo/betg;

               %Calculate 2nd order terms for metal nucleus
               if(spinm==0.0) 
                   ak = 0.0;
                   r4r5=0.0;
                   r1r2r3 = 0.0;
               else    
                    r12 = aza^2*aperp2/ak2;
                     r22 = axya2/aperp2;
               %Calculate parts of r32 and r52 which are the same
                    rp2 =(anon2*gx^2*gy^2*snphi2*csphi2)/(ak2*gsq*gperp2);
                    r32 = (rp2*aza^2*gz^2*csthe2)/(aperp2*gperp2);
                    adif = aperp2-aza^2;
                    r42=st2ct2*gz^2*gperp2*adif*adif/(ak2*gsq*gsq);
                    r52 = rp2*snthe2;
                    r1r2r3 = (r12+r22+r32)/(4*betg2);
                    r4r5 = (r42+r52)/(2*betg2);
               end

            %     start of loop over metal isotopes

            for(in=1: niso) 
                if(in>1) 
                    fwgt = tprob*wgt*pop/nline;
                     ara=arat;
                    ara2 = ara*ara;
                else    
                    fwgt= tprob*wgt/nline;
                    ara = 1.0;
                    ara2 = 1.0;
                end
           % start of loop over metal spin states

           for(isp=1: nline)
               spin=spinm-isp+1;
               spin2 = spin^2;
                hest = ho-ara*ak *spin/betg;
                hres = hest-spin2*ara2*r4r5/hest-ara2*r1r2r3*(sp2-spin2)/hest;
               if(isaveh==1)
                    if(in==iiso && isp==ilin) 
                        thedeg = the/0.01745;
                        akout = 1.0e4 * ak;
                     %Figure out later   write(15,248) thedeg, hres, akout    248: format(f6.1,2x,f7.1, 2x,f7.2)
                    end
                    if(nnit>0 || npro> 0) 
                        hyper(hint,nnit,npro) 
                    end
                    if(ishape~=7) 
                    %     this is for all cases except c60- 
                        for(int=1:nhy)
                            if(nnit>0) 
                                hhy = hres+(nnit-int+1)*akb/betg;
                            else
                                hhy = hres+(0.5*npro-int+1)*akb/betg;
                            end
                            hwgt = fwgt*hint(int)/d2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            iy1 = (hhy-out*d-hstart)/dhg + 1;
                            iyf = iy1 +2*out*d/dhg +1;
                            if(iy1<=npoint && iyf>=1) 
                                if(iy1<1) iy1 = 1; end
                                if(iyf>npoint) iyf = npoint; end
                                hoff = hstart - hhy;

                                %gaussian 1st deriv, d is peak-to-peak width
                                if(ishape==1) 
                                    fac=3.19*hwgt;
                                    for(ih=iy1:iyf) 
                                        delhd=(hoff+(ih-1)*dhg)/d
                                        y(ih)=y(ih) + fac*delhd*exp(-2*(delhd^2)); 

                                    end
                                %    lorentzian 1st deriv, d is peak-to-peak width 
                                elseif(ishape==2)
                                        fac=-.551*hwgt;
                                        for(ih=iy1:iyf)
                                            delhd=(hoff+(ih-1)*dhg)/d;
                                            y(ih)=y(ih)+fac*delhd/((0.75+delhd^2)^2); %What does ** mean???
                                        end

                               % Mixed gaussian and lorentzian first deriv
                               elseif(ishape==3)
                                   facl = -0.551*hwgt*fl;
                                    facg = -3.19*hwgt*(1.0-fl);
                                    for(ih=iy1: iyf)
                                        delhd = (hoff+(ih-1)*dhg)/d;
                                        y(ih) = y(ih)+facl*delhd/((0.75+delhd*delhd)^2);  %Again, what does ** mean???
                                    end

                               %Gaussian absorption, d is full width at half height     
                               elseif(ishape==4)
                                    fac = 0.940*hwgt*d;
                                    for(ih=iy1:iyf)
                                        delhd=hoff+(ih-1)*dhg/d;
                                        y(ih)=y(ih)+fac*exp(-2.773*delhd^2);
                                    end

                               %     lorentzian absorption, d is full width at half height
                               elseif (ishape==5)
                                    fac = 0.636*hwgt*d;
                                    for(ih=iy1:iyf)
                                       delhd = (hoff+(ih-1)*dhg)/d;
                                       y(ih) = y(ih) + fac/(1.0+4.0*(delhd^2));   
                                    end   
                               %mixed gaussian and lorentzian absorption      
                               elseif (ishape==6)
                                    facl=0.636*hwgt*d*fl;
                                    facg=0.940*hwgt*d*(1.0-fl) ;
                                    for(ih=iy1:iyf)
                                        delhd = (hoff+(ih-1)*dhg)/d ; 
                                        y(ih) = y(ih) + facl/(1.0+4.0*(delhd*delhd)) + facg*exp(-2.773*delhd*delhd);
                                    end
                                end 
                            end
                        end
                    else
                        %Do calcs for c60- distribution of g values
                        %calculate the width of the distribution at this orientation
                        %calculate 9 steps across this distribution
                        if(bwidemax<0.1) bwidemax = 0.1; end
                        bwidemin = 0.1*bwidemax

                        %bwide = bwidemin + snthe2*0.9*bwidemax
                        %the is theta in radians 
                        if(the<0.785) 
                            bwide = bwidemin;
                        else
                            bfactor = sin(2*(the-0.785))^2; %**
                            bwide=bwidemin+bfactor*0.9*bwidemax;
                        end
                        %      write(*,*) 'the, bfactor, bwide', the, bfactor, bwide 
                        stepg = bwide/3.0;
                        for(ig=1:9)
                                bg = hres+(ig-5)*stepg;
                                hoff = hstart-bg; 
                                bdif = bg-bz; 
                            if(bdif<bdifmxhf) 
                                fl=0.0;
                            else
                                fl= ((bdif-bdifmxhf)/bdifmxhf)^2; %** 
                            end 
                            d=sqrt(dz^2*((bg-bmax)/bdifmax)^2+dx^2*((bg-bz)/bdifmax)^2); %**
                            hwgt = (fwgt/(d*d))*exp(-2.773*(((ig-5)*stepg)/bwide)^2);
                            iy1 = (bg-6.0*d-hstart)/dhg + 1;
                            iyf = iy1+12*d/dhg + 1;
                            if(iy1<=npoint && iyf>=1) 
                                if(iy1<1) iy1 = 1;           end        
                                if(iyf>npoint) iyf = npoint; end
                                facl= -0.551*hwgt*fl
                                facg= -3.19*hwgt*(1.0-fl) 
                                for(ih=iy1: iyf)
                                    delhd = (hoff+(ih-1)*dhg)/d;
                                    y(ih) = y(ih)+facl*delhd/((0.75+delhd*delhd)^2) + facg*delhd*exp(-2*(delhd*delhd)); 
                                end    
                            end
                        end
                    end
               end
           end
            end
            end
        end

            %   pause 'waiting'
            if(nspex==2 && runtime==1) then

            %  store parameters for first species into variables ending with 'f'
            gxf=gx;
            gyf=gy;
            gzf=gz;
            dxf=dx;
            dyf=dy;
            dzf=dz;
            if(ishape==3 || ishape==6) then 
                flxf = flx;
                flyf = fly;
                flzf = flz;
            end
            spinmf=spinm;
            axasf=axas;
            ayasf=ayas;
            azasf=azas;
            nisof=niso;
            wgtf=wgt;
            popf=pop;
            aratf=arat;
            nnitf=nnit;
            nprof=npro;
            axbsf=axbs;
            aybsf=aybs;
            azbsf=azbs;
            theuf = theu;
            phiuf = phiu;

            %  set parameters for calculation equal to those for the second species

           gx=gxt;
           gy=gyt;
           gz=gzt;
           dx=dxt;
           dy=dyt;
           dz=dzt;

           if(ishape==3 || ishape==6) 
                flx = flxt;
                fly = flyt;
                flz = flzt; 
           end 
           spinm=spinmt;
           axas=axat;
           ayas=ayat;
           azas=azat;
           niso=nisot;
           if(niso==0) niso = 1; end
           wgt=wgtt;
           pop=popt;
           arat=aratt;
           nnit=nnitt;
           npro=nprot;
           axbs=axbt;
           aybs=aybt;
           azbs=azbt;
           theu = theut;
           phiu = phiut;

           runtime=0;
           end
    end
    
    if (nspex==2)  %  restore parameters to values for first species
        gx=gxf;
        gy=gyf;
        gz=gzf;
        dx=dxf;
        dy=dyf;
        dz=dzf;
        if(ishape==3 || ishape==6)  
            flx = flxf;
            fly = flyf;
            flz = flzf;
        end
        spinm=spinmf;
        axas=axasf;
        ayas=ayasf;
        azas=azasf;
        niso=nisof;
        wgt=wgtf;
        pop=popf;
        arat=aratf;
        nnit=nnitf;
        npro=nprof;
        axbs=axbsf;
        aybs=aybsf;
        azbs=azbsf;
        theu = theuf;
        phiu = phiuf;
    end
    
%%      if(isaveh==1) close(unit=15) ;
%       call eprplt(y,exper,npoint,legend,nlines,hstart,dhg,3)
%               %call settextwindow(r1,c1,r2,c2) 
%          301  write(*,*) ' enter 1 to change parameters '
%               write(*,*) '       2 to save current parameters '
%               write(*,*) '       3 to save simulation'
%               write(*,*) '       4 to stop '
% 
%          903  read(*,*,err=301) inext
%          902  continue
%               if (inext.lt.1 .or. inext.gt.4) then
%             write(*,*) ' value must be between 1 and 4, please re-enter '
%             go to 903
%               endif
% 
%               go to (45,501,521,999) inext
%          501  call paraout
%               go to 301
%          521  call savesim(y,npoint,hstart,hend,freq)
%               go to 301
%          999  continue
%               obox=setexitqq(qwin$exitnopersist) 
%               end
        
end

